# -*-coding: utf-8 -*-
"""
作者：TianyiFan
日期：2023-年11-月09日
用途：bender算法求解原模型
"""
import math
from gurobipy import *
from common_fuc import *
from submodel import submodel_gurobi
import para as params
import time
import numpy as np


class bender_solution:
    def __init__(self,lower_bound_value,optimal_solution):
        self.obj = 99999999
        self.solving_time = 0
        self.UB = 999999
        self.LB = 0
        self.gamma_ik = np.zeros([params.V + 2, params.K1], dtype=float)  # 主问题和子问题交互的一阶段决策变量
        self.SP_dual = []
        self.initial_LB = lower_bound_value
        self.optimal_gamma_ik = optimal_solution.gamma_ik

    def operate_benders(self):
        '''
        benders 算法的主流程
        :return:
        '''
        time_start = time.time()
        iter_index = 0
        result_change = pd.DataFrame()
        result_LB = []
        result_UB = []
        result_time = []
        result_gap = []
        result_cut = []
        while True:
            time_111 = time.time()
            # 运行MP模型，固定赋值self.gamma_ik，更新self.LB
            print()
            print(f"-------第{iter_index}次迭代----求解MP_model  结果如下----------------")
            self.build_mp_model(iter_index)
            # 将gamma_ik带入SPmodel求解，并计算本次迭代的cut_obj，并将对应的对偶变量值添加到SP_dual中,更新self.UB
            print(f"-------第{iter_index}次迭代----求解SP_model  结果如下----------------")
            if iter_index ==0:
                input_gamma_ik = self.optimal_gamma_ik
            else:
                input_gamma_ik = self.gamma_ik
            self.build_sp_model(iter_index, input_gamma_ik)
            # print(f"-------第{iter_index}次迭代----求解DSP_model 结果如下----------------")
            self.dual_SP(input_gamma_ik)
            time_222 = time.time()
            # 判断上界和下界的gap值
            gap = round((self.UB - self.LB) / self.UB * 100, 3)
            result_LB.append(self.LB)
            result_UB.append(self.UB)
            result_time.append(round(time_222-time_111,2))
            result_gap.append(gap)
            result_cut.append(self.SP_dual[iter_index].cut_obj)
            print(f"\n------第{iter_index}次迭代----benders算法概要信息----------------")
            if gap<0.01:
                print("第 {} 次迭代：下界={:.2f}  上界={:.2f} Gap={:.2f}%  耗时:{:.2f}s 停止迭代".format(iter_index, self.LB, self.UB, gap, time_222 - time_111))
                break
            else:
                print("第 {} 次迭代：下界={:.2f}  上界={:.2f} Gap={:.2f}%  耗时:{:.2f}s 继续迭代".format(iter_index, self.LB, self.UB, gap, time_222 - time_111))
                iter_index += 1

            stop_cmd = input("代码是否继续运行")
        time_end = time.time()
        self.solving_time = time_end-time_start
        self.obj = self.UB
        print("benders 分解算法的计算用时：{:.2f}s 迭代次数:{}  最终求解结果={:.2f}  迭代过程如下：  ".format(time_end-time_start,iter_index,self.obj))
        result_change['LB'] = result_LB
        result_change['UB'] = result_UB
        result_change['Gap'] = result_gap
        result_change['time'] = result_time
        result_change['Cut'] = result_cut
        result_change['delta'] = [x-y for x, y in zip(result_UB,result_cut)]
        print_pd(result_change)
        IP_model = submodel_gurobi(self.gamma_ik)
        IP_model.operate_model()
        print("将解代入IP submodel后，最优解={:.2f}".format(IP_model.obj))

    def cal_gamma_ik(self,gamma_ik):
        '''
        单独根据一阶段解计算目标函数的函数模块
        :param gamma_ik:
        :return:
        '''
        cost_ik = 0
        for i in range(params.V+2):
            for k in range(params.K1):
                cost_ik += gamma_ik[i,k]
        return cost_ik

    def build_mp_model(self,iter_idex):
        '''
        benders分解的主model
        :return:
        '''
        time_start = time.time()
        model = Model("MP_model")
        model.setParam("OutputFlag", 0)
        # model.setParam("Timelimit", 600)
        # region 定义变量
        gamma_ik = {}
        for i in range(params.V + 2):
            for k in range(params.K1):
                name5 = "gamma_{0}_{1}".format(i, k)
                gamma_ik[i, k] = model.addVar(0, 1, vtype=GRB.BINARY, name=name5)
        cut_obj = model.addVar(vtype=GRB.CONTINUOUS, name="cut_obj")
        # endregion

        # region 定义约束
        # 约束15：benders cut 的形式
        # model.addConstr(cut_obj >= self.initial_LB, name="LB_cut")
        if iter_idex>0:
            for n in range(iter_idex):
                mid_dual_set = self.SP_dual[n]
                cost1 = 0
                cost2 = 0
                # 纯参数部分
                for i in range(1,params.V+1):
                    cost1+= params.c_ui[i]*(mid_dual_set.pai_1_i[i]+mid_dual_set.pai_3_i[i])+params.c_li[i]*(mid_dual_set.pai_2_i[i] + mid_dual_set.pai_4_i[i])
                cost3 = params.C * sum(mid_dual_set.pai_8_t)
                cost = cost1+cost2+cost3
                print("第{}次迭代求解MP_model  第 {} 条cut 的 cut_cost = {:.2f}".format(iter_idex,n,cost))
                # 变量部分
                expr = LinExpr(0)
                for i in range(1,params.V+1):
                    for k in range(params.K1):
                        expr.addTerms(params.z*(mid_dual_set.pai_5_ik[i,k]+mid_dual_set.pai_6_ik[i,k]+mid_dual_set.pai_7_ik[i,k]),gamma_ik[i,k])
                model.addConstr(cut_obj>= cost+expr,name="bender_cut_{}".format(n))


        # 约束9：每层的subblock需要小于对应的params.r_ih
        for i in range(1, params.V + 1):
            expr1 = quicksum(gamma_ik[i, k] for k in range(params.K0))
            expr2 = quicksum(gamma_ik[i, k] for k in range(params.K0, params.K1))
            model.addConstr(expr1 <= params.r_ih[i][0], name="Constrain9_0_{0}".format(i))
            model.addConstr(expr2 <= params.r_ih[i][1], name="Constrain9_1_{0}".format(i))
        # 约束10：每个subblock只能分配给一个堆场
        for k in range(params.K1):
            expr = quicksum(gamma_ik[i, k] for i in range(1,params.V+1))
            model.addConstr(expr <= 1, name="Constrain10_{0}".format(k))
        # 约束11：每个time step，同一层的同一个block内，分配subblock数量不能超过3
        for t in range(params.T):
            if params.total_handle[t] >= 1:
                for g in range(params.G):
                    expr = LinExpr(0)
                    for i in range(1, params.V + 1):
                        for k in params.block_set[g]:
                            expr.addTerms(params.l_it[i][t], gamma_ik[i, k])
                    model.addConstr(expr <= 3, name="Constrain11_{0}_{1}".format(t, g))
        # 约束12：写法1--从每个subblock看，每个params.neighbor_set中，只能分配给一艘船
        for t in range(params.T):
            if params.total_handle[t] >= 1:
                for n in range(params.K1):
                    for k in range(1, len(params.neighbor_set[n])):  # 不包括自己
                        midlist = [n]
                        midlist.append(params.neighbor_set[n][k])
                        expr = LinExpr(0)
                        index_neighbor = n
                        for i in range(1, params.V + 1):
                            for index in midlist:  # --形成一个pair 对 等价于 for k in neighbor_pair
                                index_neighbor = index
                                expr.addTerms(params.l_it[i][t], gamma_ik[i, index])
                        model.addConstr(expr <= 1, name="Constrain12_{0}_({1}-{2})".format(t, n, index_neighbor))
        # 约束13：补充约束，给某艘船预留堆场的数量必定会大于总寄存量/Z 向上取整
        for i in range(1, params.V + 1):
            expr = quicksum(gamma_ik[i, k] for k in range(params.K1))
            model.addConstr(expr >= math.ceil((params.c_ui[i] + params.c_li[i]) / params.z), name="constrain13_{}".format(i))
            # print("轮船-{} 最少的subblock 与预留数量 = {}".format(i,math.ceil((params.c_ui[i] + params.c_li[i]) / params.z)))
        # 约束14：每一层堆场预留的数量数量必然会大于min(1/2预留堆场,r_ih)
        for i in range(1, params.V + 1):
            expr1 = quicksum(gamma_ik[i, k] for k in range(params.K0))
            expr2 = quicksum(gamma_ik[i, k] for k in range(params.K0, params.K1))
            model.addConstr(expr1 >= 1, name="constrain14_{}_1".format(i))
            model.addConstr(expr2 >= 1, name="constrain15_{}_1".format(i))
        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        obj.add(quicksum(gamma_ik[i, k] * min(params.t_lki[k,i],params.t_uki[k,i]) for i in range(params.V + 1) for k in range(params.K1)))
        obj.add(cut_obj)
        model.setObjective(obj, GRB.MINIMIZE)
        # endregion

        model.update()
        model.optimize()
        time_end = time.time()
        print("MP_model 求解，耗时：{:.2f}s".format(time_end - time_start))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        if model.getAttr(GRB.attr.Status) in [2, 9]:
            num_subblock = 0
            for i in range(1,params.V+1):
                for k in range(params.K1):
                    if gamma_ik[i,k].x>0.5:
                        num_subblock+= gamma_ik[i,k].x
                        self.gamma_ik[i,k] =1
                    else:
                        self.gamma_ik[i, k] = 0
            print("MP_model 每艘船预留堆场情况-------------")
            for i in range(1, params.V + 1):
                print("vessel-{}  预留subblock = ".format(i), end=" ")
                for k in range(params.K1):
                    if gamma_ik[i, k].x > 0.5:
                        print(f"{k:2} ", end="")
                print()
            # todo 检验在当前解下，加入的cut_obj值分别为多少
            if iter_idex > 0:
                for n in range(iter_idex):
                    mid_dual_set = self.SP_dual[n]
                    cost1 = 0
                    # 纯参数部分
                    for i in range(1, params.V + 1):
                        cost1 += params.c_ui[i] * (mid_dual_set.pai_1_i[i] + mid_dual_set.pai_3_i[i]) + params.c_li[i] * (mid_dual_set.pai_2_i[i] + mid_dual_set.pai_4_i[i])
                    cost3 = params.C * sum(mid_dual_set.pai_8_t)
                    cost = cost1 + cost3
                    # 带变量部分
                    cost2 = 0
                    for i in range(1, params.V + 1):
                        for k in range(params.K1):
                            cost2 += params.z * (mid_dual_set.pai_5_ik[i, k] + mid_dual_set.pai_6_ik[i, k] + mid_dual_set.pai_7_ik[i, k])*gamma_ik[i, k].x
                    totalcost=cost2+cost
                    print("第 {} 次迭代求解MP_model 加入的第{}条 cut  cut_obj = {}  参数部分 = {} 带变量部分= {:.2f}".format(iter_idex,n, totalcost,cost,cost2))
                    # print("cost_1 = {}  cost_3={}".format(cost1,cost3))
            self.LB = model.ObjVal
            # self.LB = model.ObjVal - num_subblock
            print("模型的最优解 = {1:.2f}  转化为小时={0:.2f} hours".format(model.ObjVal / 3600, model.ObjVal))
            print("预留堆场数量 = {}  原模型结果={}  LB = {}  cut_obj = {}".format(num_subblock,model.objVal-num_subblock,self.LB,cut_obj.x))
            print("预留堆场数量 = {}  原模型结果={}  LB = {}  cut_obj = {}".format(num_subblock,model.objVal-num_subblock,self.LB,cut_obj.x))

        elif model.getAttr(GRB.attr.Status) == 5:
            print("无界解")
        else:
            print("无解")
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)

    def build_sp_model(self,iter_id,gamma_ik,mode=1):
        '''
        benders分解的SP_model
        :return:
        '''
        if mode == 1:
            print("SP_model 每艘船预留堆场情况-------------")
            for i in range(1, params.V + 1):
                print("vessel-{}  预留subblock = ".format(i), end=" ")
                for k in range(params.K1):
                    if gamma_ik[i, k] > 0.5:
                        print(f"{k:2} ", end="")
                print()

        time_start = time.time()
        model = Model("SP_model")
        model.setParam("OutputFlag", 0)
        # model.setParam("Timelimit", 600)
        dual_variables = solution_dual_value()

        # region 定义决策变量
        alpha_uih = {}
        alpha_lih = {}
        beta_lik = {}
        beta_uik = {}
        for i in range(params.V + 2):
            for h in range(params.H):
                name1 = "alpha_uij_{0:>2}_{1}_".format(i, h)
                name2 = "alpha_lij_{0:>2}_{1}_".format(i, h)
                alpha_uih[i, h] = model.addVar(lb=0,vtype=GRB.CONTINUOUS, name=name1)
                alpha_lih[i, h] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=name2)
            for k in range(params.K1):
                name3 = "beta_l{0}_{1}".format(i, k)
                name4 = "beta_u{0}_{1}".format(i, k)
                beta_lik[i, k] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=name3)
                beta_uik[i, k] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=name4)
        deta_u_i = {}
        deta_l_i = {}
        abs_ui1 = {}
        abs_li1 = {}
        for i in range(1, params.V + 1):
            name6 = "deta_u_{0}".format(i)
            name7 = "deta_l_{0}".format(i)
            name8 = "abs_u_{0}".format(i)
            name9 = "abs_l_{0}".format(i)
            deta_u_i[i] = model.addVar(lb=0,vtype=GRB.CONTINUOUS, name=name6)
            deta_l_i[i] = model.addVar(lb=0,vtype=GRB.CONTINUOUS, name=name7)
            abs_ui1[i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=name8)
            abs_li1[i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=name9)
        # endregion

        # region 定义约束条件
        # 线性化约束1：最大值线性化
        pai_ih_9 = {}
        pai_ih_10 = {}
        for i in range(1, params.V + 1):
            for h in range(params.H):
                pai_ih_9[i,h] = model.addConstr(deta_l_i[i] >= alpha_lih[i, h], name="pai_ih_9_{0}_{1}".format(i, h))
                pai_ih_10[i,h] = model.addConstr(deta_u_i[i] >= alpha_uih[i, h], name="pai_ih_10_{0}_{1}".format(i, h))
        pai_i_11 = {}
        pai_i_12 = {}
        pai_i_13 = {}
        pai_i_14 = {}
        # 线性化约束2：绝对值线性化
        for i in range(1, params.V + 1):
            expr3 = quicksum(beta_lik[i, k] for k in range(params.K0))
            expr4 = quicksum(beta_uik[i, k] for k in range(params.K0))
            pai_i_11[i] = model.addConstr(abs_li1[i] >= alpha_lih[i, 0] - expr3, name="pai_i_11_{0}".format(i))
            pai_i_12[i] = model.addConstr(abs_li1[i] >= expr3 - alpha_lih[i, 0], name="pai_i_12_{0}".format(i))
            pai_i_13[i] = model.addConstr(abs_ui1[i] >= alpha_uih[i, 0] - expr4, name="pai_i_13_{0}".format(i))
            pai_i_14[i] = model.addConstr(abs_ui1[i] >= expr4 - alpha_uih[i, 0], name="pai_i_14_{0}".format(i))
        # 边界约束--可以先不要：
        # Constrain_alpha_u = {}
        # Constrain_alpha_l = {}
        # for i in range(params.V + 2):
        #     for h in range(params.H):
        #         Constrain_alpha_u[i, h] = model.addConstr(alpha_uih[i, h] <= params.c_ui[i], name="bound_alpha_u_{}_{}".format(i, h))
        #         Constrain_alpha_l[i, h] = model.addConstr(alpha_lih[i, h] <= params.c_li[i], name="bound_alpha_l_{}_{}".format(i, h))

        # 约束2/4：每艘船的unloading量守恒、loading量守恒
        pai_i_1 = {}
        pai_i_2 = {}
        for i in range(1, params.V + 1):
            expr1 = quicksum(alpha_uih[i, h] for h in range(params.H))
            expr2 = quicksum(alpha_lih[i, h] for h in range(params.H))
            pai_i_1[i] = model.addConstr(expr1 == params.c_ui[i], name="Constrain2_u_{0}".format(i))
            pai_i_2[i] = model.addConstr(expr2 == params.c_li[i], name="Constrain4_l_{0}".format(i))
        # 约束 3/5：每艘船unloading/loading到堆场的数量守恒
        pai_i_3 = {}
        pai_i_4 = {}
        for i in range(1, params.V + 1):
            expr3 = quicksum(beta_uik[i, k] for k in range(params.K1))
            expr4 = quicksum(beta_lik[i, k] for k in range(params.K1))
            pai_i_3[i]= model.addConstr(expr3 == params.c_ui[i], name="Constrain_3_u_{0}".format(i))
            pai_i_4[i]= model.addConstr(expr4 == params.c_li[i], name="Constrain_5_l_{0}".format(i))
        # 约束6：unloading 和 loading 数量在每个堆场的上界
        pai_ik_5 ={}
        for i in range(1, params.V + 1):
            for k in range(params.K1):
                pai_ik_5[i,k] = model.addConstr(beta_lik[i, k] + beta_uik[i, k] <= params.z*gamma_ik[i,k], name="Constrain6_{0}_{1}".format(i, k))
        # 约束7-8：每个堆场的loading、unloading 量上界
        pai_ik_6 = {}
        pai_ik_7 = {}
        for i in range(1, params.V + 1):
            for k in range(params.K1):
                pai_ik_6[i,k] = model.addConstr(beta_lik[i, k] <= params.z * gamma_ik[i, k], name="Constrain7_l_{0}_{1}".format(i, k))
                pai_ik_7[i,k] = model.addConstr(beta_uik[i, k] <= params.z * gamma_ik[i, k], name="Constrain8_u_{0}_{1}".format(i, k))
        # 约束13：每个time step，异层转运的数量不能超过一个阈值
        pai_t_8 = {}
        for t in range(params.T):
            expr = LinExpr(0)
            for i in range(1, params.V + 1):
                if params.l_it[i][t] == 1:
                    expr.addTerms(1 / params.handle_i[i], abs_ui1[i])
                    expr.addTerms(1 / params.handle_i[i], abs_li1[i])
            pai_t_8[t] = model.addConstr(expr <= params.C, name="Constrain13_{0}".format(t))
        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        for i in range(1, params.V + 1):
            expr1 = LinExpr(0)
            expr1.addTerms(params.time_p, deta_l_i[i])
            expr1.addTerms(params.time_p, deta_u_i[i])
            expr1.addTerms(params.time_u, abs_li1[i])
            expr1.addTerms(params.time_u, abs_ui1[i])
            expr2 = LinExpr(0)
            for k in range(params.K1):
                expr2.addTerms(params.t_uki[k][i], beta_uik[i, k])
                expr2.addTerms(params.t_lki[k][i], beta_lik[i, k])
            obj.add(expr1)
            obj.add(expr2)
        model.setObjective(obj, GRB.MINIMIZE)
        # endregion

        model.update()
        model.optimize()

        time_end = time.time()
        print("SP_model 耗时：{:.2f}s".format(time_end - time_start))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        if model.getAttr(GRB.attr.Status) in [2, 9]:
            # 赋值对偶变量
            for i in range(1,params.V+1):
                dual_variables.pai_1_i[i] = pai_i_1[i].Pi
                dual_variables.pai_2_i[i] = pai_i_2[i].Pi
                dual_variables.pai_3_i[i] = pai_i_3[i].Pi
                dual_variables.pai_4_i[i] = pai_i_4[i].Pi
                # print("p1_{}={}  p2_{}={} p3_{}={} p4_{}={}".format(i,pai_i_1[i].Pi,i,pai_i_2[i].Pi,i,pai_i_3[i].Pi,i,pai_i_4[i].Pi))
                for k in range(params.K1):
                    if  pai_ik_5[i, k].Pi!=0 or pai_ik_6[i, k].Pi!=0 or pai_ik_7[i, k].Pi!=0:
                        print("p5_{}_{}={}  p6_{}_{}={}  p6_{}_{}={}".format(i, k, pai_ik_5[i, k].Pi, i, k, pai_ik_6[i, k].Pi, i, k, pai_ik_7[i, k].Pi))
                # stop_cmd = input("代码是否继续运行")
                for k in range(params.K1):
                    dual_variables.pai_5_ik[i, k] = pai_ik_5[i, k].Pi
                    dual_variables.pai_6_ik[i, k] = pai_ik_6[i, k].Pi
                    dual_variables.pai_7_ik[i, k] = pai_ik_7[i, k].Pi

            for t in range(params.T):
                dual_variables.pai_8_t[t] = pai_t_8[t].Pi
                # print(" t = {}  pai_8_t ={}".format(t,pai_t_8[t].Pi))

            # 更新cut中的cut_obj
            dual_variables.cal_cut_obj(gamma_ik)
            self.UB = model.ObjVal
            print("第 {} 次迭代 SP_model = {:.2f}   cut_obj = {:.2f}".format(iter_id,model.ObjVal, dual_variables.cut_obj))
            self.SP_dual.append(dual_variables)

        elif model.getAttr(GRB.attr.Status) == 3:
            print("无解")
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)
        elif model.getAttr(GRB.attr.Status) == 5:
            print("无界解")

    def dual_SP(self,gamma_ik):
        time_start = time.time()
        model = Model("DSP_model")
        model.setParam("OutputFlag", 0)
        # region 定义决策变量
        pai_i_1 = {}
        pai_i_2 = {}
        pai_i_3 = {}
        pai_i_4 = {}
        pai_ik_5 = {}
        pai_ik_6 = {}
        pai_ik_7 = {}
        pai_t_8 = {}
        pai_ih_9 = {}
        pai_ih_10 = {}
        pai_i_11 = {}
        pai_i_12 = {}
        pai_i_13 = {}
        pai_i_14 = {}
        for i in range(params.V + 2):
            pai_i_1[i] = model.addVar(vtype=GRB.CONTINUOUS, name="pai_1_{}".format(i))
            pai_i_2[i] = model.addVar(vtype=GRB.CONTINUOUS, name="pai_2_{}".format(i))
            pai_i_3[i] = model.addVar(vtype=GRB.CONTINUOUS, name="pai_3_{}".format(i))
            pai_i_4[i] = model.addVar(vtype=GRB.CONTINUOUS, name="pai_4_{}".format(i))
            pai_i_11[i] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_11_{}".format(i))
            pai_i_12[i] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_12_{}".format(i))
            pai_i_13[i] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_13_{}".format(i))
            pai_i_14[i] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_14_{}".format(i))
            for k in range(params.K1):
                pai_ik_5[i, k] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_5_{}_{}".format(i, k))
                pai_ik_6[i, k] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_6_{}_{}".format(i, k))
                pai_ik_7[i, k] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_7_{}_{}".format(i, k))
            for h in range(params.H):
                pai_ih_9[i, h] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_9_{}_{}".format(i, h))
                pai_ih_10[i, h] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_10_{}_{}".format(i, h))
        for t in range(params.T):
            pai_t_8[t] = model.addVar(ub=0, vtype=GRB.CONTINUOUS, name="pai_8_{}".format(t))
        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        for i in range(1, params.V + 1):
            obj.addTerms(params.c_ui[i], pai_i_1[i])
            obj.addTerms(params.c_ui[i], pai_i_3[i])
            obj.addTerms(params.c_li[i], pai_i_2[i])
            obj.addTerms(params.c_li[i], pai_i_4[i])
            obj.add(quicksum(params.z * gamma_ik[i, k] * (pai_ik_5[i, k] + pai_ik_6[i, k] + pai_ik_7[i, k]) for k in range(params.K1)))
        obj.add(quicksum(params.C * pai_t_8[t] for t in range(params.T)))
        model.setObjective(obj, GRB.MAXIMIZE)
        # endregion

        # region 定义约束条件
        for i in range(1, params.V + 1):
            model.addConstr(pai_i_1[i] + pai_ih_9[i, 1] <= 0, name="dual_constrain_1_{}_{}".format(i, 1))
            model.addConstr(pai_i_2[i] + pai_ih_10[i, 1] <= 0, name="dual_constrain_2_{}_{}".format(i, 1))
            model.addConstr(-pai_i_1[i] - pai_ih_9[i, 0] + pai_i_12[i] - pai_i_11[i] >= 0, name="dual_constrain_3_{}".format(i))
            model.addConstr(-pai_i_2[i] - pai_ih_10[i, 0] + pai_i_14[i] - pai_i_13[i] >= 0, name="dual_constrain_4_{}".format(i))
        for i in range(1, params.V + 1):
            for k in range(params.K0):
                model.addConstr(params.t_lki[k, i] - pai_i_4[i] - pai_ik_5[i, k] - pai_ik_7[i, k] + pai_i_13[i] - pai_i_14[i] >= 0, name="dual_constrain_5_{}_{}".format(i, k))
                model.addConstr(params.t_uki[k, i] - pai_i_3[i] - pai_ik_5[i, k] - pai_ik_6[i, k] + pai_i_11[i] - pai_i_12[i] >= 0, name="dual_constrain_6_{}_{}".format(i, k))
            for k in range(params.K0, params.K1):
                model.addConstr(params.t_lki[k, i] - pai_i_4[i] - pai_ik_5[i, k] - pai_ik_7[i, k] >= 0, name="dual_constrain_7_{}_{}".format(i, k))
                model.addConstr(params.t_uki[k, i] - pai_i_3[i] - pai_ik_5[i, k] - pai_ik_6[i, k] >= 0, name="dual_constrain_8_{}_{}".format(i, k))
        for i in range(1, params.V + 1):
            expr1 = quicksum(pai_ih_10[i, h] for h in range(params.H))
            expr2 = quicksum(pai_ih_9[i, h] for h in range(params.H))
            model.addConstr(params.time_p + expr1 >= 0, name="dual_constrain_9_{}".format(i))
            model.addConstr(params.time_p + expr2 >= 0, name="dual_constrain_10_{}".format(i))
            expr3 = LinExpr()
            for t in range(params.T):
                if params.l_it[i, t] == 1:
                    expr3.addTerms(1 / params.handle_i[i], pai_t_8[t])
            model.addConstr(params.time_u + pai_i_13[i] + pai_i_14[i] - expr3 >= 0, name="dual_constrain_11_{}".format(i))
            model.addConstr(params.time_u + pai_i_11[i] + pai_i_12[i] - expr3 >= 0, name="dual_constrain_12_{}".format(i))
        # endregion
        model.update()
        model.optimize()
        time_end = time.time()
        print("DSP_model 耗时：{:.2f}s".format(time_end - time_start))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        if model.getAttr(GRB.attr.Status) in [2, 9]:
            print("DSP - obj = {} ".format(model.ObjVal))
            for i in range(1,params.V+1):
                # print("p1_{}={}  p2_{}={} p3_{}={} p4_{}={}".format(i, pai_i_1[i].X, i, pai_i_2[i].X, i, pai_i_3[i].X, i, pai_i_4[i].X))
                for k in range(params.K1):
                    if pai_ik_5[i, k].X != 0 or pai_ik_6[i, k].X != 0 or pai_ik_7[i, k].X != 0:
                        print("p5_{}_{}={}  p6_{}_{}={}  p6_{}_{}={}".format(i,k,pai_ik_5[i,k].X,i,k,pai_ik_6[i,k].X,i,k,pai_ik_7[i,k].X))
            for t in range(params.T):
                print("pai8_{} = {}".format(t,pai_t_8[t].X))

        else:
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)

class solution_dual_value:
    def __init__(self):
        self.pai_1_i = np.zeros(params.V+2,dtype=float)
        self.pai_2_i = np.zeros(params.V + 2, dtype=float)
        self.pai_3_i = np.zeros(params.V + 2, dtype=float)
        self.pai_4_i = np.zeros(params.V + 2, dtype=float)
        self.pai_5_ik = np.zeros([params.V+2,params.K1],dtype=float)
        self.pai_6_ik = np.zeros([params.V + 2, params.K1], dtype=float)  # 和一阶段决策变量有关
        self.pai_7_ik = np.zeros([params.V + 2, params.K1], dtype=float)  # 和一阶段决策变量有关
        self.pai_8_t = np.zeros(params.T, dtype=float)
        self.cut_obj = 0

    def cal_cut_obj(self,gamma_ik):
        '''
        给定一个一阶段变量，计算对偶目标函数值是否和SP的目标函数一致
        :param gamma_ik: 一阶段决策变量
        :return:
        '''
        cost1 = 0
        cost2 = 0
        for i in range(1,params.V + 1):
            cost1+= params.c_ui[i]*(self.pai_1_i[i]+self.pai_3_i[i])+params.c_li[i]*(self.pai_2_i[i]+self.pai_4_i[i])
            for k in range(params.K1):
                cost2+= params.z *(self.pai_5_ik[i,k] + self.pai_6_ik[i,k]+self.pai_7_ik[i,k])*gamma_ik[i,k]
        cost3=params.C*sum(self.pai_8_t)
        self.cut_obj = cost3+cost2+cost1
        print("cut_obj={}  cost1={}  cost2={} cost3={}".format(self.cut_obj,cost1,cost2,cost3))







